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ABSTRACT 


It Is shown that the numerical technique of Russell's momentum approach 
can be derived by using Hamilton's principle and Vance's numerical scheme. 
It results In a set of first order difference equations for solving the 
angular velocities. The method Is simple and easily programmed. The numer- 
ical examples show that the method Is also reliable. 

The algorithm Is modified next to perform the analysis of N-body sys- 
tems with closed-loop topology. To Increase the formulation flexibility, 
the equations of motion are represented by using Cartesian coordinates and 
Lagrange multipliers. The algorithm consists of two parts, Vance's scheme 
and an unconstrained minimization. The Vance's scheme Is used to find the 
angular velocities, and the unconstrained minimization Is applied to provide 
the correct angular displacements. 

The proposed scheme Is further extended to find the design sensitivity 
of an N-body system with closed-loop configuration, and to carry out the 
design optimization as well. The mmerlcal example of a small-scaled me- 
chanical system Is presented to verify the proposed formulation. Some 
aspects of future study are discussed to enhance the capability of the pro- 
posed scheme. 
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FORMULATION AND APPLICATION OF RUSSELL'S METHOD 


By 

Jean Win Hou* 

1. INTRODUCTION 

The major thrust of the Russell's method [1] for the dynamic analysis 
of multibody Is twofold. Firstly, Russell constructed a set of first order 
differential equations, uncoupled In terms of primed angular momentum. 
Secondly, the constraint forces duo to joints are eliminated In his formu- 
lation. In general, the number of first order differential equations needed 
to be Integrated are less than the number of bodies. After Integration, one 
Is left with a set of simultaneous equations for solving the angular veloc- 
ities. Russell [2] recommended the SOR (Successive Overtaxation Iter- 
ation) scheme as a solver for angular velocities. 

The Russell's method will be reformulated by using the Hamilton's 
principle and the rule of Lagrange multipliers In this report. To derive 
the Hamilton's equation for a constrained dynamic system, the variations of 
generalized coordinates and generalized velocities are treated Independently 
and the constraints are Introduced Into the derivation through the rule of 
Lagrange multipliers. The Lagrange multipliers can be Identified as con- 
straint forces. Note that the constraints for the revolute and spherical 
joints are holonomlc. Then, In section 3, the Russell's formulations for 
the N-body systems with open-loop topology are derived, and while deriving, 
the Lagrange multipliers associated with constraints are eliminated. In 
order to facilitate the development of a computer code, the equations are 

Assistant Professor, Department of Mechanical Engineering and Mechanics, 
Old Dominion University, Norfolk, Virginia 23508. 



given In the matrix and vector forms. Numerical examples of a double and 
triple pendulum are presented In section 4 to verify the aforementioned 
algorithm. 

For the N-body system with close-loop topology, the Lagrange multi- 
pliers are no longer easily eliminated. Nevertheless, the concept of 
Vance's scheme along with an unconstrained minimization scheme provides a 
very simple algorithm which Is capable of not only performing the analysis, 
but also carrying out the optimal design of such a system. The detailed 
formulation and application of this algorithm can be found In Appendix B. 

2. HAMILTON'S EQUATIONS 

For a N-degree holonomlc system, the classical approach Is to derive 
Lagrange's equations 


d 9L dL 
dt 3d 1 


0 , 


1-1....N 


from the Hamilton's principle 


( 1 ) 


ti 

6 / L dt • 0 (2) 

to 

where the Lagranglan function L Is equal to the sum of kinetic energy T 
and external work W, l.e., L ■ T + W, and 6 represents a contempora- 
neous variation. During the derivation, however. It Is assumed that opera- 

d 

tlons 6 and _ are exchangeable, l.e.. 


2 



( 3 ) 


-i«q 4 ■ «q 4 , 1"1,....,N 

dt 1 1 

In other words, the virtual velocity Is obtained by talcing the time 
derivative of the virtual displacement. Therefore, Eqs. 2 and 3 show that 

/ L(q, q, t)dt Is stationary In the family of configurations satisfying the 
N differential equations 


dq^ 

"dt 


«1 


(4) 


To make q's and q's vary Independently, the rule of multipliers asserts 
that 


ti . N dq 4 , 

/ U(q. q. t) ♦ ! X. ( — 1 - q )] dt 

t, 0 1 1 dt 1 

Is stationary for arbitrary variations of q's and q's, the X's being 
certain functions of t which are to be determined. The necessary condi- 
tions for a stationary value are given by the 2N equations 


*1 . at t 31 
dt * a^ 


1 *1, . . . , N. 


(5) 


Note that q's and the time are fixed at t 0 and t I( but not the q's* It 
can be readily verified that X's correspond to the generalized momentum 
defined In the Hamilton's principle. 

The same procedure can be extended to obtain Hamilton's equations asso- 
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dated with constrained dynamic systems. Supposing a dynamic system Is 
consistent with the following holonomlc constraints 


f^.t) - 0, 


1-1 N 

£-l,...»M, M<N 


( 6 ) 


Then, according to the rule of multipliers, there should exist £ functions 
of a(t) such that the functional 


/ ll [L(q 1f 
to 


V t) ♦ 


(^1 

dt 


q<) 


N 

+ I « f 
1 


(q< » t)] dt 


Is stationary for arbitrary variations of q's and q's. The above condi- 
tions yield 


and 


dx . ». M if. 

+ E «, JL, 1-1 N 

dt 9q 1 1 »q 1 


(7) 


X 


1 



1-1 N 


( 8 ) 


The 2N unknowns q's and q's «s well as M functions o's are to 
be determined by solving Eqs. 6-8 together. 

As to the system with nonholonomlc constraints, 


f t (q^. t) - 0 


1-1,. ...N 
£-1 H M<N 


( 9 ) 


4 


the kinematically admissible variation dq^ has' to satisfy the following 

t 

equalities. 


— i dq^ ■ 0, i«l,...,M (10) 


where d denotes a contemporaneous variation. It Is evident that Eq. 2 Is 
no longer true for the nonholonomlc system [3]. Instead one has to use the 
following equation 


/ tl dl dt - 0 (11) 

to 

Considering Eqs. 3 and 10 as the variations of constraints, the Farkas* 
lemma [4] ensures that there should exist functions A's and o's such 
that 


*i N ddq. M N Jf A 

/ [d L + E A. ( — - - dq.) ♦ I a (I -idq.)] dt - 0 
to 1 1 dt 1 1 l 3 * qi 


for arbitrary variations of q's and q's* It follows from the above condi- 
tion 




3q 1 1 



1-1 N 


(12) 


and 



( 13 ) 


Xj ■ i 

1 

Using the Eqs. 12 end 13 In conjunction with Eq. 9, the 2N unknowns, q's 
end q's, es well es M function a’s can be determined. 

From the stendpolnt of computet lonel efficiency end accuracy, formu- 
lations like Eqs. 7-8 or 12-13, for the equations of motion, ere desirable. 
Note that only N first order differential equations appear In the above 
formulations. Thus, not only the number of differential equations remains 
as N, but the potential source of error of nuierlcal Integration Is also 
minimized. However, only limited mmerlcal schemes associated with Eqs. 7-8 
or 12-13 had been developed to solve the equations of motion for dynamic 
systems. J. M. Vance [5, 6] replaced Eqs. 4 by a finite difference form and 
derived a set of finite difference equations to solve Eqs. S. 

Regarding the constrained dynamic systems, however, very few publi- 
cations are available. Numerical difficulties arise In solving q's and 

q's satisfying the constraints and In determining the corresponding multi- 
pliers. 

There are several techniques available currently for solving the 
equations of motion consistent with constraints, such as, mmerlcal 

stabilization [7, 8] and coordinate partition [9]. These methods start with 

- • 
the second order derivatives f £ (q,t) ■ 0, or first order derivative f £ 

(q, q, t) ■ 0, so that q can be a variable. In other words, an extra 

second order differential equation Is generated from each constraint. Of 

course, these approaches do not convey the original Intention of using 

Hamilton's principle which consists of first order differential equations 

only. To avoid the above difficulties, a simple way would be to eliminate 
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tht multipliers from the formulation. This Is axactly what Russell did In 
his work [1] for a dynamic systma with open-loop topology. 

3. HAMILTON'S EQUATIONS AND RUSSELL'S METHOD 
3.1 N-Body System with Cluster Configuration 
A N-body system with cluster configuration Is shown In Fig. 1. The 
X j-Y j-Z j coordinate system denotes the Inertial frame. The body-fixed 
coordinate system, X^-Yj-Zj for body 1 located at Its center. of mass. 

For the sake of easy programming, a skew-symmetric matrix I 
associated with a vector a Is defined as 



such that a x t> ■ ib * where e Is a unit vector parallel to a x S. 

It Is easy to prove the following Identities: 

? • - I (i«) 

lb ■ - (15) 

(ab)-ab-ba (16) 

Furthermore, If A Is a transformation matrix such that i • A »' and w 
Is a vector of angular velocities. It can be shown that 
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\ ■ w A 


(17) 


• • A •’ A T (18) 

,-\ bodyi*1 
i I 

body i 


Figure 1. N-8ody Cluster Geometry 
Generalized Coordinates 

Tlie position vectors of the mass centers of bodies and angular dls- 
placements of the body-fixed coordinates are taken as the generalized coor- 
dinates. The components of generalized coordinates are In terms of the 
Inertial frame. 

Generalized Forces 

As shown In Fig. 2, a force acts on the body 1 at the point P 
with the position vector P^: 



8 


*i “A, ** 


where ai* Is a vector measured from the center of mass to the point P. 
Thus, 

£l * il ♦ *1 

* «Li + h *1 

* ii ■ Ii 



Figure 2. Generalized Forces 


Replacing the differential notation d by a small increment a, one 
has 


9 




where Is the Inertia tensor of body 1 with respect to the Inertial 

C T C 

frame. The can be found equal to A^A^ where Is defined as the 
Inertia tensor with respect to the body-fixed coordinate systen and Is 
the transformation matrix between the body-fixed coordinate system and the 
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f 


Inertial frame. Thus, the variation of the total kinetic energy Is simply 
written as 


N-l _ 

«T • E (M i 4,)' 


(J^) 1 ew^ 


( 20 ) 


Constraints 

There are two sets of constraints needed to derive the Hamilton's 
equations for the constrained dynamic system. 

The first one Is a holonomlc constraint which defines the nature of 
connection between the bodies. In the following derivation, only a revolute 
joint (2 dimensional) and a spherical joint (3 dimensional) are considered. 
According to Fig. 1, the spherical (or revolute) joint provides a holonomlc 
constraint which can be expressed mathematically as 

*i ■ £L| * Sq “ .I* ■ * 0, 1*1,..., N-l (21) 

Note that all the components of the vectors In the last equation are re- 
ferred to the Inertial frame. 

The time derivative of yields 

1 ■ 0 ■ i - So - “o Vai 

d a 

After replacing the differential notation — by Incremental notation — 

dt at 

and by multiplying each term with at, one has 


11 


»S| - »So * 5, «o “ ♦ 1, 2, « ■ 0 


Thus, the variation of the constraint 6*^ >0 will yield the relation 


6 ll - % + -1 6 !o + Li 6 il “ 0 (22) 


The second set of constraints corresponds to the relations between 
displacement and velocity variations. 


d«q. 

— 1-4^ 1-0,..., N-l (23) 


and 

dee. 

. — — ■ 4w, 1-0,..., N-l (24) 

dt 


Again note that 60 ^ represents Infinitesimal angular displacement which Is 
a vector. 

For the variations of the systems, Eqs. 19, 20, 22, 23, and 24, Farkas* 

Lemma asserts that there should exist x^ (1-1 N-l), and B i 

(1-0,. ...N-l) such that for arbitrary variations 4(^, 4<[j, 46^ and 4w^: 

0 - / 1 f N f 1 [(M 1 4^ + (O*) 1 4*^] + V[fJ 4^ + (a^) 1 <6^ 

to 0 0 

+ V ii (6 ii • % * li 4 ®o + ii 4 !i) * Nfl («ai) - 4 ii] 

1 0 dt 

+ [ dT (a ®^ *' ’ dt 
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Analogous to the Hamilton's principle* £ and £ are assumed to be fixed 

at time t « t 0 and t » tj. That means ■ 0 at t * t 0 and 
t ■ tj. After Integrating by parts and collecting the corresponding terms, 
the following Identities can be obtained 


M q • a , 

1«0,...,N-1 

(25) 

1 H “1 


J w - B. , 

1*0, . . . ,N-1 

(26) 


where o and 0 can be Identified as linear momentum and angular momen- 
“1 “1 

tun (with respect to the mass center of body 1), respectively. Furthermore, 
one has 


• 

a ' 

- F 

- X 

. * °» 

1-1, ,N-1 

(27) 

"i 

"1 


1 



• 

8 ■ 

** 

• & 

F 

+ t X 

* 0, 1«1,...,N-1 

(28) 

"1 

“1 

“1 

“i"1 





• 

a 

- F + 

N f 1 X * 0 

(29) 



“0 

“0 

1 



ft 

-a F 

+ V fil X *0 

(30) 



“0 

nan 

) l “1 "1 


To reduce the nunber of degrees of freedom, two steps are necessary: 

(1) eliminating the constraint forces x , and (2) substituting the con- 

“1 

stralnts of Eq. 21, Into the formulation explicitly. 

Center of Mass 

The relations for the mass center of the whole system are defined here. 
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It can be shown that 



"f 1 H, 4t • M R, 

(31) 


N f 1 • V o 4 ■ M ft, 

0 w o - 1 

(32) 

and 

"t 1 M lll * “HR* N E 1 F, 

0 vm o “ 0 

(33) 


where M and £ denote the total mass and the position vector for the mass 
center of the whole system. The second equality In Eq. 32 follows the 
definition of as given In Eq. 25 and the last equation Is derived by 
adding Eqs. 27 and 29. Further , 

L\ a * £» 1*0,..., N-l 

or ^ “ R, + £ >i » 1*0,. ...N-l 

where r^ Is the vector between the mass centers of body 1 and the whole 
system. Eq. 21 can then be written as 

* lo • 1\ ” ■ °» 1*1,..., N-l (34) 

Then multiplying Eq. 34 with M^, summing over all the 1 and In view of 
N-l 

the fact that E M, r. * 0, It can be shown that 
0 1 

r^- - I N r 1 (dj + tj). (35) 

Ml ” 
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After substituting for j^, Eq. 34 can be rewritten as 

r i ■ - ! V M i (d^ ♦ V,) + dj, +£ 1# 1-1, . . . ,N-1. (36) 

M 1 

Furthermore, taking the time derivatives of Eq. 35 and 36, one gets 

lO " ' Y M l Go -1 + £l ^ 

and r,.-i v H, <* d, ♦ ii ) + 5o ii + ii £ 1 » ( 3Q ) 

MX 

1»1,...,N-1 

Note that £0 and rj are functions of angular velocities only. 

Constraint Forces Elimination 

The constraint forces can be eliminated from Eq. 28, by substi- 
tuting Eq. 27, 


Li Li " Gli * £j ) 1 -1, . . . N-l 

Combining the definition of and r^, as well as Eq. 33, one has 

Li “ M i £i 

M H 

* R ^ M,| r j 

N N-l 

• J. z F 4 ♦ Mi r 4 , 1-1, ...N-l (40) 

M 0 -1 1 -1 

•• 

Then Eq. 39 can be rewritten In terms of r^ as, 
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1»1,...N-1 


ii *li <"i ii> • Gi ♦ »i) £1 - li ( 5 It) 

^ • 

Since ? 4 (Mi r.) ■ . — * — - - n, t. r . , the above equation becomes 

ii <"l ii) 

M| N-l 

■ -j ' -r * <il + ii)ii ' ( Q -l) + M 1 £ i U* (41) 

S-1 ^ * 1 — 1 h o 

1“1,...,N-1 

where ■ w^ ■ w^ - t ^ w^, and is a linear combination of 

angular velocities as given by Eq. 36. 

The term 8^ + ^ Is called the primed angular momentum by 

Russell. It Is worthwhile mentioning that the right hand side of Eq. 41 Is 
a linear combination of angular velocities and the left hand side Is a qua- 
dratic function of angular velocities. 

Since only N-l equations are available In Eq. 41, one more equation has 
to be established for angular velocity wq of the central body. Making use 
of Eq. 40, the constraint force can be obtained by 


ii ■ - ii * u 




N-l 


1 ( E F. ) + M. r, + F. , 1»1,...,N-1 

MO” 1 1-1-1 


Substituting the above equality Into Eq. 30, one obtains 
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•o • ioM) 


N-l _ r - M 1 >1 4 . 

£o * ‘ ll EH, r, +_i ( I £.,) - £,] 


or 


or 


, N-l _ 

io - r ! ll < M 1 *> 


"3o4r 


* li ii) 


N-l 

*0 * \ M 1 d 1 ^1 


N-l M< N-l 
+ t 3,[- 1(1 F.) - F ] 
1 1 M 0 1 1 

-SoIo + 


N-l M, N-l 

♦ Z J [-1 ( £ FJ - FJ 

1 1 M 0 1 1 


(42) 


The above equation uses the primed angular momentum of the central body to 
calculate the angular velocity Wq. This approach 1$ different from 
Russell's formulation In that the angular momentun of the whole system Is 
established to calculate Wq. 

3.2 N-Body System with Tree Configuration 
An N-body system with tree configuration may be considered as many legs 
appended to the central body as shown <n Fig. 3. Note that a leg represents 
a structure with open-loop topology and there Is only a point joint connect* 
Ing a pair of bodies. 

So far, the equations of motion are derived only for an N-body system 
with cluster configuration. However, the derivation of equations of motion, 
discussed In the previous section, can be extended and applied to the N-body 
system with tree configuration. The derivation procedure 1$ Illustrated 
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best by deriving the equations of motion for a simple example as shown In 
Fig. 4. It Is shown that a leg which consists of four bodies Is connected 
to the central body through a single joint. In addition to the constraints 
given In Eq. 21, the constraints associated with the generalized coordinate 
of each of the legs can be expressed In terms of generalized coordinates of 
the central body as given below. By Investigating the kinematic relations 
Indicated In Fig. 4, one has 

• C|q • * 0, 1"l,...,N-2 (43) 

and 

- Sq - 4^1 - hj ■ °. j-N-l,...,N+2 (44) 

where 

4l«l * foi • 

£n * *32 * f *3 * ii3 ' idi * for 

Vi “ *13 * in * for 

and 

£( 4+2 * *34 ’ f43 + fl3 " f31 + foi 

In Eqs. 43-44, ^ and ^ are the position vectors of body 1 and the 
central body with respect to the Inertial frame. 

Note that £ Is a function of angular displacements. Thus, the varia- 
tions of h can be written as follows: 
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lx» 



body 1+1 



Figure 4. An Example of N-Body Systems with Tree Configuration. 
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a Vi " * ioi 

4 *n ■ ‘ hz a !« + ?23 " li3^ 4 1 h-i " 2oi a i-r 

a Vi "-*« aa *i ♦ <*31 " ioi> “n-i* 
and 

*Vl ’ ’ *34 ‘if 3 * <*43 ' *13> *Vl 

* ( in • ioi> *2 h-i 

Mlth the application of Hamilton's principle ind Parka's lemma, the equa- 
tions of motion for the system shown in Flq. 4 can then be derived as 

M lil*®1» 1-0,..., N+2 (45) 

J 1 *1 * £ 1 » 1-0 ** 2 (46) 

Furthermore, one also has 



2-1 " — 1 “ — 1 “ °» 


1“0,.. 

. ,N+2 

(47) 


• ** 
6 j * jS 

.1 U * h ii 

• 0, 

1-0,.. 

.,N-2 

(48) 

Vi ■ 

** — 

^*-1 £h- 1 * £oi iN- 

■1 + ioi hi 

*<*31 

* !oi J im 


- 

( hi ■ 

?01* in+2 “ 

0, 



(49) 



fj 

A + ^32 ih ' 

■ 0. 


(50) 

A# 

£n+i £#*1 

■ *i23 

_ in* hi * 

ii3 ih+i * 

<i*3 - 

il3* ih*2 “ 0| 

, (51) 
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(52) 


and 


ifH-2 


“ 1h+2 £n*2 


**34 ^2 


■ 0, 



♦ 


N*2 

r 

1 


ii 


■ o. 


(53) 


io 


•* r 

2o fo 


N-2 

r 

l 


% 


ii *3 n-i 


N*2 

( r x i ) 

N-l 


0. 


(54) 


Comparing Eqs. 45-54 with Eqs. 25-30, they are essentially Identical 
except that there are extra Lagrange multipliers Involved In Eqs. 45-54. 

Center of Hass 

The center of mass of the N-body system with tree configuration can be 
derived In the same manner as the one done In the previous section for the 
N-body system with cluster configuration. It Is not difficult to show that 

n+2 m n+2 

I #, * HR • I F. (55) 

0 - 1 “ O '" 1 


where M and JR, as defined before, denote the total mass and the position 
vector for the mass center of the entire system. By Introducing the vector 
for the relative distance between the mass center of body 1 and the 
entire system and considering the constraints, Eqs. 43-44, the relation 
between the displacements of translation and rotation can be obtained for 
each body as 
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N~2 lt+2 

Io * " jjf * M i * ii) * ^ M i (1 n-i + 


N*2 N^2 

ii • * ~ f ' *i <it ♦ »») ♦ ,* "i <4*-i ♦ ii» 




1*l,...N-2, 


N*2 Hf 2 

u ■ -i f * "i <ii ’ii) ♦ ^ "i Un-i ♦!()] 


+ in-i * -i» 


1-N-1,...,N*2 


Furthermore, taking the time derivative of the above equations, one gets 


2 u+2 

io'-jtj H, (go d, *5, £,) *t i M, (5,1, *4,)! 


ii ■ io 4 5o ^ * ii ir 


Li • io * Vi * V 


1*n-l,...N+2 


where It Is known that ^ Is a given function of angular displacements and 
velocities. As an example, ^ is obtained as 


ini * Ln+ 2 £34 * LH+1 (£43 ’ £i3) ' Sn+ 1 (£31 " lOl) 

Constraint Force Elimination 

Note that, among Eqs. 47-54, there are N+2 Lagrange multipliers In 
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total which can be completely replaced by £'s using Eq. 47. Furthermore*, 

M 

those £'s can be expressed In terms of r's* because It can be shown by 
using Eqs. 40 and 55* that 

M N+2 

a . ■ -1 I F, 1-1*. ...N+2 (62) 

1 M 0 

In addition, It Is Indicated In Eqs. 58-60 that jr 1 Is a function of angu- 
lar displacements and derivatives of angular displacements. Hence, both are 
ready to be calculated In accordance with Eqs. 47 and 62, as long as the 
values of £, e^, and are valid. The position of mass center of the 
entire system Jl, can be Integrated from Eq. 55. As to the angular dis- 
placements and their derivatives, they also can be Integrated based on the 
equations associated with primed angular momentums. Since the derivations 
of those equations are similar to the ones presented In the last section, 
oniy the primed angular momentun of the central body is derived here as an 
Illustration. 

Substituting Eq. 47 Into Eq. 54 for the Lagrange multipliers, one could 
obtain 


N-2 N+2 N-2 

lo ‘ 2<£o ‘ \ -i -i " ?N+1 U ) + \ ^ 1\ 

mm 

Furthermore, the term can be replaced by jr^ 
Eq. 62, 


N+2 # 

+ < 3 kl . ( E Oi ) * 0 
N-i N-l - 1 


with the application of 
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* ~ N-2 _ _ H+2 N-2 Mi N+2 

lo • JoEo • 5 li£t - Vi £t> * 5 % (jf { f|> 


N+2 M 1 N+2 


N-2 


N+2 


* Vi ^ r ( 5 *.»♦ I M, u Vi M < Vi r, 


- 0 


Using the equality, uv ■ uv + uv, the above equation can be rearranged to 

# 

provide an equality for the primed angular momentum of the central body as: 


N-2 


N+2 


So + \ M 1 -£l + M 1 Vl-1 

N-2 N+2 

• lo£o + \ llli + Vl ^ -i* 

N-2 M, N+2 N+2 M. N+2 

- E (J. E FJ - Vi [ I — ( r £.)] 

1 1 M 0 1 N-l M 0 ^ 

N-2 


+ e Mi (tj^ -5^) r, 


N+2 


+ n *1 M 1 0M-i3d * 5o3n-1> ±4 


(63) 


where the left side of the equality denotes the total derivative of the 
primed angular momentum of the central body. Note that the preceding equa- 
tion Is a perfect form suitable for applying Vance's nunerlcal scheme [5, 
6]. The numerical Implementation of such a scheme for the dynamic systems 
of concern Is to be discussed In the next section. 

4. NUMERICAL ALGORITHM ANO EXAMPLES 
The numerical Implementation of Russell's formulation Is discussed 
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hereafter. Numerical examples are also presented to verify the derived 
formulation's and numerical algorithm. 

For simplicity, Eqs. 33, 38, 41-42 and Eq. 63 can be represented sym- 
bolically by 


fc (L» £» H) * £ (£• £» £) ( 64 ) 

where angular velocity jv ■ £ and ]> denotes position vector such as 4, 
<1 or £ defined with respect to the Inertia frame. It Is evident that £ 
depends on £. Furthermore, the primed angular momentum £ Is a linear 
function of jv, and J5, on the other hand. Is a nonlinear function of j*. 
To Implement Eq. 64 numerically, the major step Is to approximate the dif- 
ferential operator by a difference operator [6]. In the numerical examples 
discussed below, the trapezoidal rule is used, l.e., 

Vl *f-t S in> + G ( !4l+l*£n*l'2n+l>I <65 > 

where the subscript denotes the time grid point. Since Eq. 65 represents a 
set of nonlinear algebraic equations for w n+ j, an Iterative scheme is 
required to find w^. The detailed numerical algorithm Is described In 
Fig. 5. 

There are two examples, a double pendulum and a triple pendulum. The 
results obtained according to the Russell's formulation are compared with 
those calculated by directly employing Lagrange's equations of motion with- 
out Introducing Lagrange multipliers. The latter can be done when the mini- 
mum set of generalized coordinates are selected to describe the system. For 
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the examples studied, the angular displacement of each pendulun Is defined 
as a generalized coordinate. The Lagrange's equations of motion constitute 
a set of simultaneous second-order equations which are generally nonlinear: 


— <— > 
dt )£ 


— ♦ I (a. i. t) 


Mote that the term — Is subjected to the total differentiation. Thus, 

a i 

the Lagrange's equations of motion can be replaced by a set of finite dif- 
ference equations which are used to solve £, following the same numerical 
Implementation as used In the Russell's formulation. The detailed formula- 
tion of Lagrange's equations of motion, as well as the Russell's formula- 
tion, for a double and a triple pendulum are give In Appendix A. 

The mass, length and the moment of inertia at the mass center of each 
pendulum are assigned as 1, 1 and 0.08333 units, respectively. A two 
unit force Is applied at the free end of the systems and is always normal to 
the pendulum. The Initial configuration of the pendulum Is straight In the 
vertical position as shown In Figures 6-7. The numerical results calculated 
with at * 0.05 seconds are listed in Tables 1 and 2 for the double and 
triple pendulum, respectively. The history of motion Is depicted in Figures 
6-7 as well. The algorithm performs well. In general. It takes two to four 
Iterations for jv and to get convergence (e * 10* 5 ) at each time grid 
point. For the examples studied. It Is Indicated that the Russell's formu- 
lation and the Lagrange's equations of motion provide essentially the same 
results. Nevertheless, the Russell's formulation requires more CPU time 
because the introduction of joint reactions (Lagrange multipliers) and the 
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Figure 7a. Initial Configuration 



a.a 2.m t.n i.a «.a ia.a u.a 

Tire GRID (SEC) ■10" 

Figure 7b. Changes of Angular Velocities 



Table 1. Numerical results for a double pendulum 
a. Results based on the Lagrange's equations of motion 
(CPU time - 1.17 sec. on DEC 10) 


Time (Sec.) 

0! (Rad.) 

O 2 (Rad.) 

(Rad/Sec.) 

#2 (Rad/Sec.) 


O.OOOOOE+OO 

0.00000E+00 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.10 

•0.8554 5E -02 

0.42832E-01 

-0.17076E+00 

0.8561 3E+00 

0.20 

-0.33358E-01 

0.16995E+00 

-0.31859E+00 

0.16755E+01 

0.30 

-0.68293E-01 

0.37174E+00 

-0.36078E+00 

0.2331 4E+01 

0.40 

-0.99210E-01 

0.62794E+00 

-0.23739E+00 

0.27714E+01 

0.50 

-0.11053E+00 

0.92259E+00 

0.26391E-01 

0.31194E+01 

0.60 

-0.89484E-01 

0.12531E+01 

0. 40967E+00 

0.34990E+01 

0.70 

-0.23541E-01 

0.16261E+01 

0.92792E+00 

0.39764E+01 

0.80 

0.1Q224E+00 

0.20531E+01 

0. 16093E+01 

0.45796E+01 

0.90 

0.30438E+00 

0.25461E+01 

0.24503E+01 

0.52892E+01 


0.59483E+00 

0.31115E+01 

0.33567E+01 

0.60147E+01 


b. Results based on the Russell's formulation 
(CPU) time « 3.19 sec. on OEC 10) 

Time (Sec.) 

Oi (Rad.) 

O 2 (Rad.) 

wj (Rad/Sec.) 

M 2 (Rad/Sec.) 

0.00 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.10 

-0.85545E-02 

0.42832E-01 

-0.17076E+00 

0.85613E+00 

0.20 

-0.33358E-01 

0.16995E+00 

-0.31859E+00 

0.16755E+01 

0.30 

-0.68293E-01 

0.37174E+00 

•0.36078E+00 

0.23314E+01 

0.40 

-0.99210E-01 

0.62794E+00 

-0.23739E+00 

0.27714E+01 

0.50 

-0.11053E+00 

0.92259E+00 

0.26391E-01 

0.31194E+01 

0.60 

-0.89484E-01 

0.12531E+01 

0.40967E+00 

0.34990E+01 

0.70 

-0.23541E-01 

0.16261E+01 

0.92792E+00 

0.39764E+01 

0.80 

0. 10224E+00 

0.20531E+01 

0.16093E+01 

0.45796E+01 

0.90 

0.30438E+00 

0.25461E+01 

0.24503E+01 

0.52892E+01 

1.00 

0.59483E+00 

0.3U15E+01 

0.33567E+01 

0.60147E+01 
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position vector of mass center of the entire system complicates the compu- 
tation. 

5. DISCUSSION ANO CONCLUSION 

The Vance's scheme has been successfully applied to analyze the dynam- 
ics of mechanical systems. The same scheme Is also extended to perform the 
design sensitivity analysis and the optimum design. The proposed algorithm 
Is simple and easily programmed. 

There Is no need for this algorithm to evaluate the higher order deriv- 
atives of equations of constraints. Moreover, If the acceleration of the 
dynamics system Is of no concern. It Is also not necessary to find the time 
derivative of the mass matrix. Thus, the proposed algorithm relieves , to 
some extent, the complexity of formulation and computation of mechanical 
system dynamics. 

It Is noted that the time Increment, at. In this algorithm plays an 
Important role, not only for the numerical error and stability but also for 
the convergence rate of the Iterative scheme, because the extrapolation of 
the current state variable with a smaller at provides a better estimate of 
the new state variable for starting the Iterative computation at the new 
time grid. 

The example of the slider-crank mechanism discussed In Appendix B shows 
the success of the application of the proposed algorithm for solving a 
small-scale problem. In order to analyze the large-scale problem, the pro- 
posed algorithm can be upgraded by replacing the unconstrained minimization 
scheme Introduced In Eqs. B.9-11 by a more efficient scheme designed for the 
large nonlinear problem. As an example, the Fletcher-Reeves's conjugate 
gradient algorithm could be a promising substitute. 
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In order to fully develop Russell's momentum approach, the following 
tasks are also proposed: 

1. Extend the derivation to Include flexibility, and to extend to 
systems with 3-0 configuration. 

2. Extend the derivation to Include complicated joint conditions, such 
as Joint friction. 

3. Investigate the application of the sparse matrix technique to 
Improve the computational efficiency of solving angular veloc- 
ities. 

In suamary, while a simple strategy Is proposed and successfully Imple- 
mented to analyze and optimally design a small-scaled mechanical system with 
dlfferentlal/algebralc equations, further study Is required to enhance the 
algorithm's rigorousness and versatility. 
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APPENDIX A 


EQUATIONS OF NOTION FOR DOUBLE ANO TRIPLE PENDULUM 



A.l. Lagrange's Equations of Motion 

Double Penduluw 

Definitions of notations are given In Figure A.l. 


A 



Figure A.l. A Double Pendulu*. 
The equations of notion for a double pendulim are 


J o“o * loio 2o * *1 loi^i * 


J £i ' "A ioii " m ii£i£o2o * -\- 


and 


where J Q and are the moments of Inertia of pendulums 0 and 1 with 
respect to points A and B, respectively. 

Triple Pendulum 

Definitions of notations are given In Figure A. 2. 



The equations of motion for a triple pendulum are 


JqSo - (nyi^) 

- (mja^ + n.^ ♦ m{^) ♦ ^F. 


and 


J lil * ("illio + m 2-A ^ ^0 ’ m 2ll^l— 1 * *2^2^2 
- « 1 i 1 l 1 * 0 !!o + ™ 2 ^1 ( £o*0 * + il-* 


JgWg - " | gj g*Q 3iQ - •2fi£A 
• (IqWq ♦ i^) + £2 F 


where Jq, ^ and Jg are the moments of inertia of pendulums defined with 
respect to the points A, B and C, respectively. 

A. 2. Russell's Equations of Motion 
A. 2.1. Double Pendulum 

Definitions of notations are given in Figure A.3. 



Figure A. 3. A Double Pendulum. 


Russell's equations of motion for a double pendulum are 


J oSo ♦ ll ("ill) 


Vo ♦ "1 <g&4&> k - vAr ££ 1* 


m^ im^ 


and 


♦ ii <Ml> 


<V*i> L * — !o l * -li 


n, l 4fl> 2 


Vl-i 


where 


and 


h * - — G ot ♦Siii* + S<& + ii£i‘ 

*n!l^ 


t n* <v m w 

ii • £ “ !!i£i * -ill - £aJ!i 


Triple Pendulua 

Definitions of notations are given in Figure A.4. 








Russell's equations of motion for a triple pendulum are 


J 2*2 + 

J A * io 2a * ?i2’Iio ) 22 

and 


IgF *^1 +lzi?2* 


JqSo 





2o<£o 


(?oo + i)i *2i " ^Do^ioi ] h 
-(loo^i* — * 


where 


and 


tO 






Vo2oo 
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21 

22 


IIH 

- - *1 loo * (— * »i5 ii 
"to 

m l m 2 • L 

■ - "to ioo * — k* 
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100 “ Vbo * *cioo ” looV 

101 “ iioi * *0*01 * loi^o* 

iio- liiio - SOio - h£v 
il2 “ -lil2 " -ill 2 ■ !lA» 

igl - w 2 t 21 - Wg£ 21 - 
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APPENDIX B 


DYNAMIC ANALYSIS AND OPTIMIZATION 
OF 

CONSTRAINED MECHANICAL SYSTEMS* 


♦Submitted to the International Journal for Numerical Methods In Engineering 
for the consideration of publication. 



B.I. INTRODUCTION 


For a complex mechanical system. It has been Indicated In the literature 
[1, 2] that the equations of motion may be presented by using Cartesian 
coordinates and lagranglan multipliers, resulting In a system of mixed 
differential and algebraic equations. This approach greatly Increases the 
formulation flexibility, because It does not rely on the engineer's Intuition 
to determine a set of Independent variables. 

Several nunerlcal methods can be found In the literature Cl, 2] to solve 
these differential -algebra 1c equations. One method Is to differentiate the 
equations of constraints and append them to the equations of motion. This 
expanded system of equations Is solved for the lagranglan multipliers and 
accelerations. The Integrations of accelerations provide good predictions for 
the velocities which must be subsequently corrected based on the equations of 
constraints. The velocities and displacements can be treated in a similar 
way. A large number of state variables Is usually encountered In this 
approach which may become a prohibitive problem to be solved. To overcome 
this difficulty, a coordinate partition scheme [3] Is Implemented at each time 
step to sort out the Independent and dependent coordinates based on the 
equations of constraints. Only the Independent coordinates are to be 
Integrated. Recently, the singular value decomposition scheme has been 
Introduced [4, 5] to select a set of composite coordinates as Independent 
coordinates which constitute a hyperplane tangent to the equations of 
constraints. The numerical study shows that the singular value decomposition 
method Is quite promising. 

An easily programmed algorithm Is presented In this paper. The algorithm 
consists of two parts: (a) a scheme Introduced by Vance [6, 71 and (b) an 
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unconstrained minimization. Vance's scheme has been limited to solve 
unconstrained equations of motion. However, In this study, Vance's scheme, 
along with an unconstrained minimization algorithm Is used to analyze the 
dynamics of a constrained mechanical system. An unconstrained minimization Is 
Implemented to correct the state variables by minimizing the constraint 
violations. The proposed scheme Is further extended to find the design 
sensitivity of constrained mechanical systems by the adjoint variable 
technique [ 8 ] and to carry out the problems of optimization as well. 

B.II. DYNAMICS ANALYSIS 

A mechanical system Is defined as a system that consists of bodies with 
Inertias and elements without Inertias such as control force, damper, etc. To 
define a mechanical system, one may assign a body-fixed coordinate for each 
body and Introduce the equations of contralnts to describe the kinematic 
relations between the bodies. Each body has either three degrees of freedom 
corresponding to a two dimensional configuration or six degrees of freedom 
corresponding to a three dimensional configuration. In this way, the kinetic 
energy and the external work of each Individual body can be easily 
established. Based on Hamilton's principle and the theorem of lagranglan 
multipliers, the equations of motion of a whole system can then be derived as 

It (-5 " - * - to* i» (B - l) 

and a set of constraints as 

• ( 4 . t) • 0 . (B.2) 

The total kinetic energy T, a quadratic form of velocities, is the sum of the 
kinetic energy of each body. The terms £, £ and ^ denote the generalized 
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coordinates, velocities and forces, respectively. In this study, the 
constraint vector. ®, Is limited to a set of holonomic constraints. Note 
that the generalized coordinates g and largranglan multipliers \ are the 
unknowns In the above equations of motion. The numerical Implementation of 
this system of differential /algebra 1c equations Is discussed next. 


For simplicity, Eqs. B.l and B.2, can be represented symbolically by 

jjT- ' (ig) T £ * £ i t) (B.3) 

where the mo men tun term 8T/&£ 1$ equal to M£. At this stage, one may 
Introduce Vance's scheme to approximate the differential operation by a finite 
difference operator. As an example. If the trapezoidal method Is employed, 
Eq. B.3 can be replaced by a set of finite dlffeience equations defined at 
time (n+1) At and nAt: 


8® T 


8® T 


1 I n+1 * 1 l n + *7 + In + h + ln+1^* 

The preceding formula Is usually a nonlinear equation with roots, and 

To find them, the simple linear Iteration Is sufficient and 
convenient. Let and be good Initial estimates of solutions 

3 i, + l and After rearrangement, Eq. B.4 may be rewritten as the 

following recursive form for jth Iteration: 


^n+l* * ^ n+1 * ^^n+1^ 


a® 


(j) T 


At 


" M ( V in + ^7 in * * ^7 * ^nll (B ’ 1 * * * 5) 

Although the above equations become linear equations with and 

^in+l» they are unable to solve both and because the 
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number of unknowns Is larger then the number of equations. Nevertheless, with 
the help of equations of constraints, one obtains the following Identity by 
differentiating Eq. B.2, 

M • » 

It may also be written In a recursive form defined at the time .equal to tp+j. 


5» (j) #(1+1) » (J) 


(B.6) 


(J+l) 


The last equation along with Eq. B.5 provides a matrix form to solve £ 
and simultaneously: 


M 

"n+1 

f n (J> 

( ^n+l 


-f n U) 

( *#n+l 


b® T 


•(M) 
i n+1 

at v (j+i) 
“7- n+1 


KtJ i + £] n +^G ( ^ 1 

bft (j) 

^"5? n+1 


(B.7) 


The leading coefficient matrix can be proved to be positive definite 

Oft 

provided that the rows of are linearly Independent [3], Thus, the 

existence and uniqueness of and ^ \ are ensured. The new 

value of can then be obtained by numerically Integrating £^j*j , 

for Instance, by using the trapezoidal method. 




(B.8) 
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• f 1+1 I <• 

Since the generalized velocities £ ^ ore not Independent, the generalized 
coordinates obtained by Integrating or* not kinematically 

permissible. In other words, may not satisfy the equations of 

constraints, l.e., 9 To consistent with the 

equations of constraints, an unconstrained minimization scheme Is proposed to 
simply reduce the deviations of 9 t), l.e., 

Min ♦ ••(£. t) T • (q, t) (B.9) 

Vl 

where the design variable Is The Initial estimate of Is 

provided by the direct Integration of Eq. B.8. 

There are many methods available to carry out the unconstrained 

minimization defined In Eq. B.9. Numerical results presented In section 4 are 

obtained by a recursive quadratic programming algorithm [9], called the 

linearization method, which has been proved to be globally convergent. More 

* 

specifically, the new value of 5 Is obtained by modifying the current value 
of £ In the following way: 


3 . 


* 2a 9 1 


a* 

4 


(B.10) 


where the parameter a Is a step size determined In such a way that the cost 

* 

<J» 0 Is always reduced for the Improved value of £ , l.e., 

*0 tf > < ♦ 0 ( £ } ' Be II £ T * ^ M 2 (B * 11) 


where e Is a given constant, usually defined ss 0.1, and the notation 
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|| • || denotes the L 2 norm. The computations of Eqs. 8.10 and 8.11 
constitute an Iterative process to be terminated whenever the value of 
|| • || becomes very small. After the value of Is updated by the 

optimum solution of the unconstrained mlnlmzatlon, Eq. 8.9, the Iterations 
between Eqs. B.7-9 continue until both g n<|>1 and reach the convergence 

criteria: 



( 8 . 12 ) 


I A ( ^ 1} 
I - n+1 


X (J) 

- n+1 


I < • 


where the notation | • | denotes the L* norm and the e Is a given small 
constant. Once the convergence Is achieved, the computation moves to the next 
step and the Iteration starts again. 

The numerical algorithm Is summarized as follows: 

Step 1: Start with Initial conditions and X » 0. 

Step 2: Select the Initial values for and X^|j. 

Step 3: Solve the matrix equation (8.7) for 1^}^ and 

Step 4: Calculate the Initial estimate of by using equation 

8 . 8 . 

Step S: Update the value of by carrying out the unconstrained 

minimization, Eqs. 8.9*11, In order to correct the constraint 
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deviation* 

Step 6: Check the convergence criteria defined In equation B.12, If the 

convergence Is achieved, move to Step 2 with n • n*l; otherwise, 
j ■ j+1, move to Step 3. 

There are some remarks worthwhile mentioning here. First, If It Is difficult 
to find the Initial conditions for all of the dependent £ and £ In Step 1, 
the unconstralnted minimization scheme given In Step 5 with Independent 
variables fixed can be used to obtain accurate dependent £ and £ . Second, 
the generalized acceleration £* can be calculated directly by rearranging Eq. 
B.3 as 

to T 

£ * H ^ [A £ ♦ (^) k ♦ £ (£, £, t) ] 

Compared to the method Introduced In references 2 and 3, the proposed approach 
avoids the complicated and time-consuming process of deriving the second order 
derivative of equations of constraints, l.e., 4* . 

B.III. DESIGN SENSITIVITY ANALYSIS 

The design sensitivity analysis of a system with differential and 
algebraic equations associated with mechanical system dynamics has recently 
been a subject of study [10]. Two approaches have been discussed In the 
literature. One Is the direct differentiation method [11]. The other Is the 
adjoint variable technique [8, 12]. 

The dlfferentlal/algebralc equations for the dynamics of a constrained 
mechanical system can be rewritten as follows: 
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and 

• (b, 3, t) ■ 0 

with the Initial conditions glvon as 

at t • 0 (B.13) 

i ■ at t ■ 0 . 

Tht kinetic ontrgy T Is glvon as a quadratic function of l.o., 

T ■<y^ T M b, t) j for a symmetric matrix M. 

It Is vory common In tht optlnal design formulation to havo the cost or 
constraint written In a functional form as 

t 

♦ • / F (b, £, £, dt 

o “ 

The task for the design sensitive analysis Is to obtain the design 
derivative of a with respect to the design variable b. The design 

variation of 4*. 64> Is derived as 



or after Integrating by parts, 

*■/’ ♦ Eli - ] a’t «** * r? a* i; <»•>*) 

0 - * 0 £ j 

where the terms with apostrophes denote the variations due to the 
perturbations of design variables, 6b. It Is revealed In Eqs. B.I and B.2 
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that the relation between design variable b and state variable £ Is highly 
nonlinear. Thus, It Is difficult to find £' explicitly In terms of ]>. 
Nevertheless, the adjoint variable technique provides an alternative In which 
£* Is not required to be defined analytically. 

To begin the adjoint variable technique, one pre-multl piles two arbitrary 
vectors p(t) and ^(t) to Cqs. B.l-2, and Integrates the products over the 
time period (0, t) to get the Identities: 

M- p T M £ - p T |I - jJq - ji T (*^) T \] dt ♦ p T M £ \l - 0 (B. 15 ) 


and 


^ T 

J v‘ * dt « 0 
o “ “ 

The design variation of Eq. B.1S Is derived as 


/ { • P T (M 2 *~ jf T (M £ 5b ” i T m i' 


•T, 


•T 




(8 .16) 





a 


-£ t [» /l l. a a-i T tw/S, , 
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(• 

* *1 


a > T A'} dt 


♦ 




♦ j* T (M 5) b 6b | - 0 

where the subscripts denote derivatives and where the terms with a bar on top 
are not subject to the differentiation. Integrating the terms with In 
the preceding equation by parts, one obtains 

Ma ,T ^<"£> ♦C(§>I A d. t *[aJ a d. t 

■ (M »£ Ji " (-3^) Ji " Ji " t£ ® •*! 

■ 6 J> t A> »b T k * ^ + £»J ** ■ 1> T »b jrf 

- x ,T dt 

" £ T M S' 1 1 * (|j) » £ I* 1 1 ‘ i/ £• ^ S' IJ 

+ “ T (m i , » ± s' i' to ■ / <m &>• 1 • »i: * ° 


Furthermore, the design variation of Eq. B.16 provides 
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/ (J ®,^ £* ♦ v* • , ^ to) dt » 0. 

Adding the last two equations to Eq. 8.14 and grouping the corresponding terms 
together, one Is able to express the design variation of the given functional 
4> symbolically as 

6 4 • / (A* 6b — hJ £* + u T ® T , A )dt + boundary terms (B.17) 

o “ -a “ 

where = (M , T b ji ♦ (|I) , T b u Q, b T \x 

* £ »b - + -*b T - + F »b T 
and 

^q s *3? (M + t(||)» T £ iil't + t2*q T l*3»t 

' < M A>> a £ ‘ £ 

• tS T 2* a 3» a ii + ®* T a * + F » a T • < F ^>» t T 

It Is noted that the design variation 6<t> Is a linear functional of 
variations, 6J>, j}' and A', and that the vectors £ and still remain 

unspecified. One may then assign values for u and v_ so that A^ * 0 and 
u ■ 0. In rt.ier words, the adjoint variables u and are defined so 
as to eliminate the Influence of unknowns £* and A* In the formulation of 
6<|). The condition A^ ■ 0, along with * 0, can be arranged and 

written as follows: 
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(M £) ♦ £ ■ R (J>, jj, J, £, £, £, t) (8.18) 

and 

*,^£•0 (8.19) 

where R symbolically denotes the rest of the terms defined In A^. It Is 
evident that the adjoint equations 8.18 and 8.19 are linear functions In terms 
of adjoint variables £ and £ and they represent a mixed system of 

dlfferentlal/algebralc equations. More specifically, Eq. 8.19 provides a set 
of linear constraints on £, and v serves as a vector of lagranglan 

multipliers. 

Following a similar procedure one can Investigate the boundary terms 
shown In Eq. 8.17 and determine the terminal conditions of the adjoint 
variable £ *1n such a way that the Influence of design variations £' and £' 
Is eliminated from the boundary terms of 64». Note that the boundary terms 
of 64> have the design variations £' and £* defined at both t * 0 and t * 
x. If the Initial conditions of all the Independent and dependent 
coordinates, as well as velocities, are given explicitly, for Instance: 

£ * % (£) , at t ■ 0 

i ■ a, (w ■ *‘ • ’ °- 

then the design variations £* and at t « 0 are found without difficulty 
as 

£* * at t * 0 
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Hence, the boundary terms of 64 can be rearranged to obtain: 

Boundary terms - 

- { [m; T M ♦ »i T (jl) » £ ♦ J* T Q» £ “ j* T < M (^) 

' ii” M (5^) " ± (M 1> *b) * 6 £l t»0 

- {[»*^ M ♦ j*^ (<jl) ± £,^ * (M 3) ,^] q_^ ♦ M £*} t 

♦£ T (M |), b 6b|J (B.20) 

On the other hand, while the boundary conditions are known only for the 
Independent coordinates and velocities l.e.. 


£l *3o at t- ° 

il *io — at t *°» 


(B.21) 


the equations of constraints defined at t*o should be used In order to find 
the. design variations of the dependent quantities. The design variations of 
Eq. B.21 and Eq. B.2 at t*> provide the following matrix equation 
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8 ® 8 ® 8 # 


(B.22) 


Which can be used to determine the design variation explicitly In terms 

of 6b at t-0 Provided that the matrix b®/bj^ Is not singular. Further- 
more, the design variations of £} and 3' 0 can be obtained In a similar way 
by constructing Eq. B.22 based on the design variation of & at t « 0. 
Therefore, the boundary terms of 64» are still able to be written In the form 
of equation B.20, although only the Independent coordinates and velocities are 
given at t^O. 

Next the boundary terms of 6* defined at t»x are Investigated. The 
best way to avoid calculating the unknown variations £ and at the 
terminal time x Is to specify the terminal conditions of adjoint variable 
li it t ■ t so that the terms associated with and £' can be dropped. To 
achieve this. It Is sufficient to obtain the following Identifies from Eq. 
B.20, 


M £ ♦ [(|I) , T | + Q T .£ - (M £) ,^] T £ - 0, at t-x 
and 

M i* ■ 0 , at t»x 

Because of the positiveness of the mass matrix M, It is simply concluded from 
the above conditions that p(x) ■ jl (x) - 0. Based on these two terminal 

conditions, along with the adjoint equations, Eqs. B.18-19, the adjoint 
variables u(t) and Wt) can then be determined uniquely In the entire period 
of time (0, x). Note that the same mmerlcal scheme used for dynamic 
analysis, Eqs. B.l-2, can be applied here to solve the adjoint equation 



B. 15 


numerically, though a bettor schema could be Implemented to take advantage of 
the linearity of the adjoint equations. 

Finally, with the knowledge of £, v, £and X, the design variation of 
the cost functional 6*, Eq. B.14, can be expressed as a linear functional of 
the perturbation of the design variable, 6b, 

6* » J T aJ 6b dt ♦ ^ 6bf 0 ♦ j* T (M £) , b 6b 1 1 

0 — 

where Is defined In Eq. B.17. While the Initial condition of all the 
Independent and dependent coordinates, as well as velocities, are given 
explicitly, the term is defined as 



Otherwise, a similar form can still be obtained based on Eq. B.22. 

B. IV. NUMERICAL EXAMPLE 

A modified slider-crank mechanism with one degree of freedom Is studied 
here as an example to validate the numerical algorithm presented previously. 
This mechanism Is composed of two linkages, the crank AB and the connecting 
rod BC, as shown in Fig. B.l. While each of the hinge joints A and B 
entertains two constraints, the joint point C Is forced to slide along the x- 
axls. With a torque H applied at the joint A, the system Is subjected to a 
planar motion. The definitions of body-fixed coordinates, as well as the 
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vectors Aj and d 1# ara Indicated In Fig. B.2. The total kinetic energy of 
the system Is obtained as 


2 

Z 

1 


1 

7 


(Mi ai 


♦ J i w*) 


The kinematic constraints for joints A. B and C are given as. 


♦i 2 ♦ £i “ 0 

♦a 5 ii + ii * ^ " 0 

♦ 3 S J T (Rj + dg) - o 


(B.23) 


where the unit vector 2 1s parallel to the Y-dlrectlon of the Inertia frame. 
The constraint 4 3 means that there Is no Y-component of joint C's movement 
at any time. 


Analysis 

Based on Hamilton's principle and theory of lagranglan mullpllers, a 

system of eleven equations can be set up In the form of Eq. B.7 for the 

slider-crank mechanism. The unknowns to be solved are the six degrees of 

freedom, Rj, J^, w 3 and w ? , as well as the five lagranglan multipliers, 

and \ 3 associated with the equations of constraints, Eq. B.23. The 

diagonal components of the 6 x 6 mass matrix M are Mj, Mj, M 2 , M 3 , Jj and 

J 2 » The fifth component, equal to the given torque H, Is the only non- zero 

element In the forcing term jj. The detailed formulation of the Jacobian, 

» 

, Is given In the Appendix. 

With the non-dimensional data: Mj»M 2"1, H»2, Jj ■ 0.08333, J 2 ■ 0.33333, 

jtjl - |dj| « 0.5, and |x 2 | • 1^1 ■ 1» the slider-crank mechanism of 
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concern Is analyzed on the DEC-10 system with double precision. The time step 
At Is set as 0.005 seconds. It takes 25.5 CPU seconds to simulate the 

motion for the time period of one second. The convergence criteria, Eq. 8.12, 
for the coordinates and velocities are set as 10~ 5 . It takes at most three 
Iteratlves at each time grid to achieve the given convergence requirement. As 
for the unconstrained minimization scheme for the coordinate correction. It 
also takes at most three Iterations at each time grid to achieve the 

convergence given as 

11*11 < u >' 6 

Some of the miner leal results are listed In Table B.l. The last column 

contains the deviation of the constraint, + 3 * regarding the sliding joint 
C. Furthermore, the results, obtained by using a commerlcally available 
program, DAOS [13], are also listed In Table B.l for the purpose of 

comparison. A good agreement between the purposed scheme and DADS Is 
observed. 

Design Sensitivity Analysis and Optimization 

The optimization of slider-crank mechanism studied here Is to find the 
control torque H(t) so that the motion of the sliding joint C can follow a 
desired path n(t). This example falls Into the category of Inverse 
dynamics. Nevertheless, an optimization formulation Is set up to approximate 
the best torque profile H(t) In terms of given functions. That Is 

T 9 

Min 4 > 0 - / [R 2 + Ag - n(t)r dt 
H(t) 0 

where the control torque Is expressed In terms of given functions Nj(t) and 
design variables bj as 
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H(t) ■ £ b i H^t). (B.24) 

The function Nj(t) can be a polynomial or trigonometric function. The 
desired path* n ■ fo x * ti y ] T , is given as 

n x (t) » a ♦ 10 (p-a) t 3 - 15 (p - a)t 4 + 6 (p-o) t S , 0 < t < 1, 
n y (t) ■ 0 

Note that the path n x (t), starting from the Initial position a to the final 
position p, entertains zero velocities and accelerations at both t«0 and 
t*T. The design variation of the cost function Is derived as 

t 8Q j 

64*0 * / Jf * 6H 
0 

• 

where the forcing term Q Is given as (0, 0* 0* 0* H* 0) T , and where 
1 

6H ■ z Nj(t) • 6qj because of equation B.24. Thus* 

t 

64» 0 ■ E[J Mg • Nj (t) dt] • 5qj 
0 

5 1 T (B.25) 

where A Is defined as a sensitivity vector. The adjoint variable jx Is the 
solution of the following adjoint equations: 

jj”. + 1 " * + 2 + + JI ) (B<26) 


and 
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q £ • 0 (B.27) 

with terminal conditions; ji(t) - £ (t) ■ 0. The detailed formulation of the 
terms In Eq. B.26 can be found In the Appendix. 

The accuracy of the design sensitivity analysis Is Investigated next. 
First, the control torque Is fixed as H»2. The design sensitivity vector 
calculated by the analytical equation Is compared with the one calculated by 
the direct finite difference method. The result plotted In Fig. B.4 shows a 
very good agreement between the aforementioned twu methods when the 
perturbation of the design variable Is up to 101. Second, the control torque 
Is described by a quadratic function H(t) * a 0 + ajt + ajt 2 with three 
coefflcents as design variables. The components of the design sensitivity 
vector provided by Eq. B.2S, as well as those calculated by the finite 
difference method are listed In Table B.2. The results are obtained based on 
a 0 ■ aj ■ 82 ■ 1 and 21 perturbation for each design variable. It Is 
Interesting to observe from the design sensitivity vector that to Increase the 
design coefficient a 0 Is more beneficial In terms of the reduction of error In 
the path generation than to Increase a^ and a 2 . 

Once the accurate design sensitivity Is produced, any gradient-based 
mathematical programming can be used to generate the optlimm solution 
Iteratively. The following nunerlcal results are obtained by using a 
recursive quadratic program called the linearization method [9], The control 
torque Is assuned to have one of the following forms: 

H( t) ■ a 0 , 

H(t) ■ a 0 t ajt, 

H(t) • a 0 ♦ ajt ♦ a 2 t 2 . 



H(t) ■ si nut ♦ b 2 cositt 

where the coefficients, a 0 , a 2 , bj and b 2 are treated as design variables. 

Corresponding to differently prescribed torque functions, the optimum 
solutions and their associated data are listed In Table B.3, and plotted In 
Fig. B.5 as well. Note that none of the trial torque functions are able to 
establish a path pattern so as to make the terminal velocity and acceleration 
approach zero. Finally, to demonstrate the stable performance of the optimum 
algorithm, the convergence progressions of the cost function and the 
convergence criteria, l.e., L 2 norm of the design gradient, are also plotted 
In Fig. B.6 only for the quadratic control torque. 
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Tabic B.l. Numerical Results of Analysis 


Time 

Prop. 

Algorithm 

DAOS 


♦3 

6j(Rad) 

"1 

OjIRad) 

"1 

0.1 

0.01499 

0.29993 

0.01697 

0.29987 

0.32xl0" 19 

0.2 

0.05992 

0.59759 

0.06346 

0.59698 

0.14xl0“ 17 

0.3 

0.13410 

0.88239 

0.13822 

0.88085 

0.13xl0* 16 

0.4 

0.23524 

1.1324 

0.23958 

1.1294 

0.45xl0* 16 

0.5 

0.35865 

1.3254 

0.36317 

1.3201 

0.38xl0' 16 

0.6 

0.49826 

1.4578 

0.50234 

1.4522 

0.29xl0* 16 

0.7 

0.64877 

1.5473 

0.65248 

1.5415 

0.14xl0’ 17 

0.8 

0.80714 

1.6189 

0.81052 

1.6136 

0.35xl0’ 16 

0.9 

0.97270 

1.6946 

0.97595 

1.6906 

0.17xl0“ 15 

1.0 

1.1468 

1.7927 

1.1501 

1.7908 

0.34xl0' 15 
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Table B.2 Design Sensitivity Analysis of the Sllder-Cranker Mechanism 

With Quadratic Control Torqua: 

H(t) • a 0 ♦ ajt ♦ a 2 t* 



id 

Adjoint Variable 

m 

Finite Difference 



Technique 

Method 

Ooslgn Gradients 

eq. 

[ttejtbaj) -e(a f ))Ma 1 






0 

•0.22357 

-0.219 

2.04* 

*0 

5ij 

-0.05845 

-0.059 

0.94S 


^o 

*r 2 


•0.02389 


-0.0245 


2.93* 



Table B.3 Optimum Control Torques of the Slider-Crank Mechanism 
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ADDENDUM TO APPENDIX B 
W 

1. The Jacobian, ^ : 



where the first subscript denotes the number of the body and the second 
subscript represents the component along the specified direction. 


2. The terms on the right side of Eq. B.26, 

(a) ( i T sf i^’q 

i ° 

o 

0 

0 

A lx ^ h * A ly h * d lx '5 *3 + d ly ^5 X 4 
' *2x ^6 *3 " A 2y ^ X 4 * d 2y ^6 X 5 


(b) 2 (Rg + dg) 1 .^ • (Kg ♦ dg - n) 

r 8. . 


2(R 

+ d 

- n ) 

2x 

2x 

X 

2(R ♦ d 

2y 2y 

) 

•o 

CM 

« 

(R 

♦ d 

2y 

2x 

2x 

2 d 

(R 

♦ d 

2x 

2y 

2y 


V 
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Fig. B.4. Design Sensitivity Analysis of the Slider-Crank Mechanism with 
Constant Control Torque: H * 2. 


3000 


B.31 


i.* 

8 

C 

3 


o> 

3 

O’ 

o 



m x l* 33NV1SI0 


o 

Lil 

<S) 

LLJ 

z 


Fig. B.5. The Optimal Paths for Different Torque Functions 




Fig. B.6. Numerical Convergence History o^ the Slider-Crank 
Mechanism ( with quadratic torque ) 






